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C/} Abstract In the present work we have selected a collection of statistical and 

^ mathematical tools useful for the exploration of multivariate data and we present 

them in a form that is meant to be particularly accessible to a classically trained 
mathematician. We give self contained and streamlined introductions to princi- 
pal component analysis, multidimensional scaling and statistical hypothesis test- 
ing. Within the presented mathematical framework we then propose a general 
exploratory methodology for the investigation of real world high dimensional 
datasets that builds on statistical and knowledge supported visualizations. We ex- 
emplify the proposed methodology by applying it to several different genomewide 
DNA-microarray datasets. The exploratory methodology should be seen as an em- 
bryo that can be expanded and developed in many directions. As an example we 
point out some recent promising advances in the theory for random matrices that, 
if further developed, potentially could provide practically useful and theoretically 
. ^ well founded estimations of information content in dimension reducing visualiza- 

tions. We hope that the present work can serve as an introduction to, and help to 
o3 stimulate more research within, the interesting and rapidly expanding field of data 

exploration. 
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1 Introduction. 

In the scientific exploration of some real world phenomena a lack of detailed 
knowledge about governing first principles makes it hard to construct well-founded 
mathematical models for describing and understanding observations. In order to 
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gain some preliminary understanding of involved mechanisms and to be able to 
make some reasonable predictions we then often have to recur to purely statistical 
models. Sometimes though, a stand alone and very general statistical approach 
fails to exploit the full exploratory potential for a given dataset. In particular a 
general statistical model a priori often does not incorporate all the accumulated 
field- specific expert knowledge that might exist concerning a dataset under con- 
sideration. In the present work we argue for the use of a set of statistical and 
knowledge supported visualizations as the backbone of the exploration of high 
dimensional multivariate datasets that are otherwise hard to model and analyze. 
The exploratory methodology we propose is generic but we exemplify it by ap- 
plying it to several different datasets coming from the field of molecular biology. 
Our choice of example application field is in principle anyhow only meant to be 
reflected in the list of references where we have consciously striven to give refer- 
ences that should be particularly useful and relevant for researchers interested in 
bioinformatics. The generic case we have in mind is that we are given a set of ob- 
servations of several different variables that presumably have some interrelations 
that we want to uncover. There exist many rich real world sources giving rise to 
interesting examples of such datasets within the fields of e.g. finance, astronomy, 
meteorology or life science and the reader should without difficulty be able to pick 
a favorite example to bear in mind. 

We will use separate but synchronized Principle Component Analysis (PCA) 
plots of both variables and samples to visualize datasets. The use of separate 
but synchronized PCA-biplots that we argue for is not standard and we claim 
that it is particularly advantageous, compared to using traditional PCA-biplots, 
when the datasets under investigation are high dimensionsional. A traditional 
PCA-biplot depicts both the variables and the samples in the same plot and if 
the dataset under consideration is high dimensional such a joint variable/sample 
plot can easily become overloaded and hard to interpret. In the present work we 
give a presentation of the linear algebra of PCA accentuating the natural inherent 
duality of the underlying singular value decomposition. In addition we point out 
how the basic algorithms easily can be adapted to produce nonlinear versions of 
PCA, so called multidimensional scaling, and we illustrate how these different 
versions of PCA can reveal relevant structure in high dimensional and complex 
real world datasets. Whether an observed structure is relevant or not will be judged 
by knowledge supported and statistical evaluations. 

Many present day datasets, coming from the application fields mentioned 
above, share the statistically challenging peculiarity that the number of measured 
variables (p) can be very large (10 4 < p < 10 10 ), while at the same time the num- 
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ber of observations (AO sometimes can be considerably smaller (10 1 <N< 10 3 ). 
In fact all our example datasets will share this so called "large p small N" char- 
acteristic and our exploratory scheme, in particular the statistical evaluation, is 
well adapted to cover also this situation. In traditional statistics one usually is 
presented with the reverse situation, i.e. "large N small p", and if one tries to 
apply traditional statistical methods to "large p small N" datasets one sometimes 
runs into difficulties. To begin with, in the applications we have in mind here, 
the underlying probability distributions are often unknown and then, if the num- 
ber of observations is relatively small, they are consequently hard to estimate. 
This makes robustness of employed statistical methods a key issue. Even in cases 
when we assume that we know the underlying probability distributions or when 
we use very robust statistical methods the "large p small N" case presents difficul- 
ties. One focus of statistical research during the last few decades has in fact been 
driven by these "large p small N" datasets and the possibility for fast implemen- 
tations of statistical and mathematical algorithms. An important example of these 
new trends in statistics is multiple hypothesis testing on a huge number of vari- 
ables. High dimensional multiple hypothesis testing has stimulated the creation 
of new statistical tools such as the replacement of the standard concept of p-value 
in hypothesis testing with the corresponding q-value connected with the notion of 
false discovery rate, see lfT2l . lfT3l . Il48l . Il49l for the seminal ideas. As a remark 
we point out that multivariate statistical analogues of classical univariate statis- 
tical tests sometimes can perform better in multiple hypothesis testing, but then 
a relatively small number of samples normally makes it necessary to first reduce 
the dimensionality of the data, for instance by using PCA, in order to be able to 
apply the multivariate tests, see e.g. lfT4l . Il33l . Il32ll for ideas in this direction. In 
the present work we give an overview and an introduction to the above mentioned 
statistical notions. 

The present work is in general meant to be one introduction to, and help to 
stimulate more research within, the field of data exploration. We also hope to 
convince the reader that statistical and knowledge supported visualization already 
is a versatile and powerful tool for the exploration of high dimensional real world 
datasets. Finally, "Knowledge supported" should here be interpreted as "any use 
of some extra information concerning a given dataset that the researcher might 
possess, have access to or gain during the exploration" when analyzing the visu- 
alization. We illustrate this knowledge supported approach by using knowledge 
based annotations coming with our example datasets. We also briefly comment 
on how to use information collected from available databases to evaluate or pre- 
select groups of significant variables, see e.g. [fTT | . [fT6l . 0TI .[[42 | for some more 
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far reaching suggestions in this direction. 



2 Singular value decomposition and Principal com- 
ponent analysis. 

Singular value decomposition (SVD) was discovered independently by several 
mathematicians towards the end of the 19'th century See [47J for an account 
of the early history of SVD. Principal component analysis (PCA) for data anal- 
ysis was then introduced by Pearson ll25l in 1901 and independently later devel- 
oped by Hotelling [|24l . The central idea in classical PCA is to use an SVD on 
the column averaged sample matrix to reduce the dimensionality in the data set 
while retaining as much variance as possible. PCA is also closely related to the 
Karhunen-Loeve expansion (KLE) of a stochastic process [|29l . 041 . The KLE of 
a given centered stochastic process is an orthonormal L 2 -expansion of the process 
with coefficients that are uncorrected random variables. PCA corresponds to the 
empirical or sample version of the KLE, i.e. when the expansion is inferred from 
samples. Noteworthy here is the Karhunen-Loeve theorem stating that if the un- 
derlying process is Gaussian, then the coefficients in the KLE will be independent 
and normally distributed. This is e.g. the basis for showing results concerning the 
optimality of KLE for filtering out Gaussian white noise. 

PCA was proposed as a method to analyze genomewide expression data by 
Alter et al. LI J an d has since then become a standard tool in the field. Super- 
vised PCA was suggested by Bair et.al. as a regression and prediction method 
for genomewide data flHl , fl9j , ifTTl . Supervised PCA is similar to normal PCA, 
the only difference being that the researcher preconditions the data by using some 
kind of external information. This external information can come from e.g. a re- 
gression analysis with respect to some response variable or from some knowledge 
based considerations. We will here give an introduction to SVD and PCA that 
focus on visualization and the notion of using separate but synchronized biplots, 
i.e. plots of both samples and variables. Biplots displaying samples and variables 
in the same usually twodimensional diagram have been used frequently in many 
different types of applications, see e.g. EDI , ETTl . E2l and lfT5l , but the use of sep- 
arate but synchronized biplots that we present is not standard. We finally describe 
the method of multidimensional scaling which builds on standard PCA, but we 
start by describing SVD for linear operators between finite dimensional euclidean 
spaces with a special focus on duality. 
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2.1 Dual singular value decomposition 

Singular value decomposition is a decomposition of a linear mapping between 
euclidean spaces. We will discuss the finite dimensional case and we consider a 
given linear mapping L : R N — > R p . 

Let ei, e2, • • • , e>j be the canonical basis in R N and let fj, f2, . . . , f p be the canon- 
ical basis in R p . We regard R^ and R p as euclidean spaces equipped with their 
respective canonical scalar products, (-,-)r# and (v)rp, in which the canonical 
bases are orthonormal. 

Let L* : R p — > R N denote the adjoint operator of L defined by 

(L(u),v) R , = (u,L*(v)V ;ueR w ;veR". (2.1) 

Observe that in applications L(et), k = 1,2, . . . ,N, normally represent the arrays of 
observed variable values for the different samples and that L*(fj), j = 1,2, .. . ,p, 
then represent the observed values of the variables. In our example data sets, the 
unique pxN matrix X representing L in the canonical bases, i.e. 

X jk = (fj,L(e fe )) R p ;j= l,2,...,/7; k= l,2,...,iV, 

contains measurements for all variables in all samples. The transposed N x p 
matrix X T contains the same information and represents the linear mapping L* in 
the canonical bases. 

The goal of a dual SVD is to find orthonormal bases in R N and R p such that 
the matrices representing the linear operators L and L* have particularly simple 
forms. 



We start by noting that directly from (2.1 ) we get the following direct sum 
decompositions into orthogonal subspaces 

R N = KerL@lmV 

(where KerL denotes the kernel of L and ImL* denotes the image of If) and 

R p =ImL®KerL\ 

We will now make a further dual decomposition of ImL and ImL* . 

Let r denote the rank of L : R w — > W, i.e. r = dim {ImL) = dim (ImL*). The 
rank of the positive and selfadjoint operator L* oL : H N — > R w is then also equal 
to r, and by the spectral theorem there exist values Ai > A2 > ■ • ■ > A,- > and 
corresponding orthonormal vectors u x ,u 2 , . . . ,u r , with u k G R^, such that 

L*oL(u k ) = \ 2 u k ;Jfc=l,2,...,r. (2.2) 
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If r < N, i.e. dim (KerL) > 0, then zero is also an eigenvalue for L* oL : — > 
R N with multiplicity N — r. 

Using the orthonormal set of eigenvectors {u 1 , u 2 , . . . , u r } for L* o L spanning 
ImL* , we define a corresponding set of dual vectors v 1 , v 2 , . . . , v r in R p by 



L(u k ) =: V 



From (2.2) it follows that 



L*(v k ) = A,u k ;k 



1,2, 



1,2,. 



and that 



;k=l,2, 



(2.3) 

(2.4) 
(2.5) 



LoL*(v k ) = \ 2 v k 

The set of vectors {v 1 , v 2 , . . . , v r } defined by (2.3 ) spans ImL and is an orthonor- 
mal set of eigenvectors for the selfadjoint operator LoL* : R p — > W. We thus 
have a completely dual setup and canonical decompositions of both R^ and W 
into direct sums of subspaces spanned by eigenvectors corresponding to the dis- 
tinct eigenvalues. We make the following definition. 

Definition 2.1 A dual singular value decomposition system for an operator pair 
(L, L*) is a system consisting of numbers X\ > X2 > ■ ■ ■ > K > and two sets 
of orthonormal vectors, {u^u 2 , . . . ,u r } and {v 1 ^ 2 , . . . ,v r } with r = rank(L) = 
rank(L*), satisfying ( 2.2\ -(2.5) above. 

The positive values X\ , A2, . . . , X r are called the singular values of(L, L*). We 
will call the vectors u , u 2 , . . . , u r principal components for ImL* and the vectors 
v 1 , v 2 , . . . , v r principal components for ImL. 

Given a dual SVD system we now complement the principal components for 
ImL*, u^u 2 , . . . ,u r , to an orthonormal basis u^u 2 , . . . ,u N in R w and the prin- 



cipal components for ImL, v 1 , v 2 , 
RP. 

In these bases we have that 



, v r , to an orthonormal basis v 1 , v 2 . . . , v p i 



(vU(u k )) R P = (r(vj),uV 

This means that in these ON-bases L : R^ - 
p x N matrix 

D 




in 



ifj ''^ r (2.6) 
otherwise . 

R p is represented by the diagonal 

(2.7) 
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where D is the rxr diagonal matrix having the singular values of (L, L*) in de- 
scending order on the diagonal. The adjoint operator L* is represented in the same 
bases by the transposed matrix, i.e. a diagonal N x p matrix. 

We translate this to operations on the corresponding matrices as follows. Let 
U denote the N x r matrix having the coordinates, in the canonical basis in R N , for 
u , u 2 , . . . u r as columns, and let V denote the pxr matrix havin g the coordinates, 



in the canonical basis in R p , for v^v 2 , . . .v r as columns. Then (2.6) is equivalent 
to 

X = VDU T and X T = UDV T . 

This is called a dual singular value decomposition for the pair of matrices (X, X T ). 

Notice that the singular values and the corresponding separate eigenspaces for 
L* o L as described above are canonical, but that the set {u 1 , u 2 , . . . , u r } (and thus 
also the connected set {v 1 , v 2 , . . . , v r }) is not canonically defined by L* o L. This 
set is only canonically defined up to actions of the appropriate orthogonal groups 
on the separate eigenspaces. 

2.2 Dual principal component analysis 

We will now discuss how to use a dual SVD system to obtain optimal approxima- 
tions of a given operator L : — > R p by operators of lower rank. If our goal is 
to visualize the data, then it is natural to measure the approximation error using a 
unitarily invariant norm, i.e. a norm || ■ || that is invariant with respect to unitary 
transformations on the variables or on the samples, i.e. 

IILII = ||VoLo[/|| for ally and £/ s.t. V*V = Id andU*U = Id . (2.8) 



Using an SVD, directly from (2.8) we conclude that such a norm is necessarily a 
symmetric function of the singular values of the operator. We will present results 
for the L 2 -norm of the singular values, but the results concerning optimal approx- 
imations are actually valid with respect to any unitarily invariant norm, see e.g. 
Il35ll and ll40l for information in this direction. We omit proofs, but all the results 
in this section are proved using SVDs for the involved operators. 

The Frobenius (or Hilbert- Schmidt) norm for an operator L : R N — > R p of 
rank r is defined by 

\\ L \\f : = > 

where X^,k= 1 , 2, . . . , r are the singular values of (L, L*). 
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Now let J'nx,, denote the set of real n x n matrices. We then define the set of 
orthogonal projections in R" of rank s < n as 



m ■- 



{TieJ^nxn ; n* = 11; non = n; rank(Tl) =s}. 



One important thing about orthogonal projections is that they never increase the 
Frobenius norm, i.e. 

Lemma 2.1 Let L:R N — >W be a given linear operator. Then 

||IIoL||f<||L||f forallTle&P 

and 



||LoII||i? < \\L\\ F for allTle 
Using this Lemma one can prove the following approximation theorems. 
Theorem 2.1 Let L : R N — >RP be a given linear operator. Then 

sup ||n p oLoIT^II^ = sup ||rioL||/r = 

TIP € &>! ; U N € ne 

{min(s,/-) ^1/2 
£ *}) (2-9) 
k=i J 



and equality is attained in (2.9) by projecting onto the min(5,r) first principal 
components for ImL and ImL*. 

Theorem 2.2 Let L:R N — y RP be a given linear operator. Then 



inf ||L-lFoLorr|| F = inf ||L-IIoL|| F = 

up € &>f ; u N e &>n ne 

{max(.s,r) ^ ^'^ 

£ Xl \ (2.10) 
fc=min(.s,r) + l J 



and equality is attained in {2.10) by projecting onto the min(s,r) first principal 
components for ImL and ImL*. 
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We loosely state these results as follows. 
Projection dictum 

Projecting onto a set of first principal components maximizes average projected 
vector length and also minimizes average projection error. 

We will briefly discuss interpretation for applications. In fact in applications 
the representation of our linear mapping L normally has a specific interpretation in 
the original canonical bases. Assume that L(e k ), k= 1,2,... ,,/V represent samples 
and that L*(fj), j = l,2,...,p represents variables. To begin with, if the samples 
are centered, i.e. 

£L(e k )=0, 

k=\ 

then \\L\\p corresponds to the statistical variance of the sample set. The basic 
projection dictum can thus be restated for sample-centered data as follows. 
Projection dictum for sample-centered data 

Projecting onto a set of first principal components maximizes the variance in the 
set of projected data points and also minimizes average projection error. 
In applications we are also interested in keeping track of the value 

X jk = (f i ,L(e k )). (2.11) 

It represents the j'th variable's value in the fc'th sample. 

Computing in a dual SVD system for (L, L*) in ( |2.11[ ) we get 

X jk = M (e k , u 1 ) (fj, v 1 ) + • ■ ■ + A,(e k , u r ) (fj, v r ) . (2. 12) 



Now using (2.3) and (2.4) we conclude that 

Xjk = f (e k ,L* (v 1 ))(f j ,L(u 1 )) + ■ ■ ■ + i (e k ,L* (v r ))(fj,L(u r )) . 
Finally this implies the fundamental biplot formula 

X jk = f^yXL^y) + ■ • ■ + l(L(e k ), v r )(L* (fj),u r ) . (2.13) 
We now introduce the following scalar product in R r 

(a,b)x := ^-aibi-\ V^-a r b r ;a,beR r . 

A\ Ay 
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Equation (2.13) thus means that if we express the sample vectors in the basis 
v 1 , v 2 , . . . , v r for ImL and the variable vectors in the basis u 1 , u 2 , . . . , u r for ImL* 
, then we get the value of Xj% simply by taking the (-, -^-scalar product in W 
between the coordinate sequence for the k'th sample and the coordinate sequence 
for the j'th variable. 

This means that if we work in a synchronized way in W with the coordinates 
for the samples (with respect to the basis v^v 2 , . . . , v r ) and with the coordinates 
for the variables (with respect to the basis u^u 2 , . . . ,u r ) then the relative posi- 
tions of the coordinate sequence for a variable and the coordinate sequence for a 



sample in W have a very precise meaning given by ( |2.13 ). 



Now let S C { 1 , 2, . . . , r} be a subset of indices and let |S| denote the number 
of elements in S. Then let n£ : R p — > W be the orthogonal projection onto the 
subspace spanned by the principal components for ImL whose indices belong to 
S. In the same way let IT^ : — > R N be the orthogonal projection onto the 
subspace spanned by the principal components for ImL* whose indices belong to 
S. 

If L(et), k = 1,2, . . . ,iV represent samples, then IT^ oL(ek), k = 1,2, , iV, 

represent S -approximative samples, and correspondingly if L*(f J ), j = 1,2, . . . ,p, 
represent variables then Il^oL*(fj), j = 1,2,...,/?, represent S-approximative 
variables. 

We will interpret the matrix element 

Xj : =(fj,n£oL(e k )) (2.14) 

as representing the j'th S-approximative variable's value in the 
fc'th S-approximative sample. 



By the biplot formula (12. 1 3b for the operator IT^ oL we actually have 



X %= I ^-(^e k ),v m )(r(fj),u m ). (2.15) 

meS A,n 

If |S| < 3 we can visualize our approximative samples and approximative variables 
working in a synchronized way in Rl 5 l with the coordinates for the approximative 
samples and with the coordinates for the approximative variables. The relative 
positions of the coordinate sequence for an approximative variable and the coor- 
dinate sequence for an approximative sample in Rl 5 l then have the very precise 



meaning given by (2. 15 ) 



Naturally the information content of a biplot visualization depends in a crucial 
way on the approximation error we make. The following result gives the basic 
error estimates. 
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Theorem 2.3 With notations as above we have the following projection error es- 
timates 



p N 

EEft*-*M 2 =El*l 2 (2-i6) 

;'=lfe=l i£S 

sup \X jk - Xf h | < sup | Xi | . (2.17) 

We will use the following statistics for measuring projection content: 

Definition 2.2 With notations as above, the L? -projection content connected with 
the subset S is by definition 



a.JS) := ^|| ! 
Li=i \M\ 



2 



We note that, in the case when we have sample centered data, 052(5) is precisely 
the quotient between the amount of variance that we have "captured" in our pro- 
jection and the total variance. In particular if 052(5) = 1 then we have captured 



all the variance. Theorem 2.3 shows that we would like to have good control of 
the distributions of eigenvalues for general covariance matrices. We will address 
this issue for random matrices below, but we already here point out that we will 
estimate projection information content, or the signal to noise ratio, in a projec- 
tion of real world data by comparing the observed L 2 -projection content and the 
L 2 -projection contents for corresponding randomized data. 



2.3 Nonlinear PCA and multidimensional scaling 

We begin our presentation of multidimensional scaling by looking at the recon- 
struction problem, i.e. how to reconstruct a dataset given only a proposed covari- 
ance or distance matrix. In the case of a covariance matrix, the basic idea is to 
try to factor a corresponding sample centered SVD or slightly rephrased by taking 
the square root of the covariance matrix. 

Once we have established a reconstruction scheme we note that we can apply 
it to any proposed "covariance" or "distance" matrix, as long as they have the 
correct structure, even if they are artificial and a priori are not constructed using 
euclidean transformations on an existing data matrix. This opens up the possibil- 
ity for using "any type" of similarity measures between samples or variables to 
construct artificial covariance or distance matrices. 
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We consider a p x N matrix X where the N columns {xi, . . . , x^} consist of 
values of measurements for N samples of p variables. We will throughout this 
section assume that p>N. We introduce the N x 1 vector 

l = [l,l,...,l] r , 

and we recall that the N x N covariance matrix of the data matrix X = [xi, . . . , x^] 
is given as 

c( Xl , . . . ,x N ) = (x - ixi i r ) r (x - Ixi l T ) . 

We will also need the (squared) distance matrix defined by 

D^(xi,...,x N ) := |xj-x k | 2 ; j,k = 1,2, . . . ,N. 

We will now consider the problem of reconstructing a data matrix X given 
only the corresponding covariance matrix or the corresponding distance matrix. 
We first note that since the covariance and the distance matrix of a data matrix X 
both are invariant under euclidean transformations in R p of the columns of X, it 
will, if at all, only be possible to reconstruct the pxN matrix X modulo euclidean 
transformations in R p of its columns. 

We next note that the matrices C and D are connected. In fact we have 

Proposition 1 Given data points xi,X2, . . . ,xn in R p and the corresponding co- 
variance and distance matrices, C and D, we have that 

D/fc = C jj + C» - 2C j k . (2. 1 8) 

Furthermore 

C Jk = ^ t (»y + D tt ) - \» jk D ta , (2. 19) 

LLy i=\ z Ziv i,m=l 



or in matnxnotatwn, 



Proof. Let 



i N 

y i :=x i --£x j ; i= 1,2,. ..,N. 

" 7=1 
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Note that 



Cjk = yjyk and D/fc = |yj-yk|' 

and that lj} =l yj = 0. 

Equality (|2.18[) above is simply the polarity condition 



Moreover, since 



D^=| yj | 2 + |y k | 2 -2C^. 



N N 

£C /jt = and £C 7l = : 



by summing over both j and k in (2.21 ) above we get 



A 7 i A 7 



7=1 



On the other hand, by summing only over j in (2.21 ) we get 



N N 

LD ; , = A/-|v k | 2 +£| yj | 

7=1 7=1 



Combining ( |2.22| ) and ( |2.23[ ) we get 

1 N 



1 w 

7=1 7,«=1 



Plugging this into (2.21 ) we finally conclude that 



D 7*=^I( D *7+D 



A? 



ik) 



;=1 



/V 2 . 



jk ■ 



',7=1 



This is (2.19). q.e.d 



(2.21) 



(2.22) 



(2.23) 



Now let ^NxN denote the set of all real N xN matrices. To reconstruct apxN 
data matrix X = [xj, . . . , xn] from a given N x N covariance or N x N distance 
matrix amounts to invert the mappings: 



4> : R p x ■•• x R p 3 (x 1 ,x 2 ,...,x N ) H> C(x 1 ,x 2 ,...,x N ) £^ NxN , 
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and 



W-.RPx ■■■xRP 3 (x!,x 2 ,...,x N ) h-*.D(xi,x 2 ,...,x N ) eJZ NxN . 

In general it is of course impossible to invert these mappings since both 4> and *F 
are far from surjectivity and injectivity. 

Concerning injectivity, it is clear that both 4> and *F are invariant under the 
euclidean group E(p) acting on R p x • ■ • x R p , i.e. under transformations 

(x 1} x 2 ,...,x N ) i y (5x! + b,5x 2 + b,...,5x n + b) , 

where b G R p and S G O(p) . 

This makes it natural to introduce the quotient manifold 

(RPx---xRP)/E(p) 

and possible to define the induced mappings <J» and *F, well defined on the equiv- 
alence classes and factoring the mappings 4> and *P by the quotient mapping. We 
will write 

<t> : ([xi,x 2 ,...,x N ]) !->■ C(xi,x 2 ,...,x N ) 

and 

W: ([xi,x 2 ,...,x N ]) h^D(x 1 ,x 2 ,...,x n ). 

We shall show below that both 4> and *F are injective. 

Concerning surjectivity of the maps 4> and or 4> and we will first de- 
scribe the image of 4>. Since the images of 4> and *F are connected through Propo- 
sition [TJ above this implicitly describes the image also for It is theoretically 
important that both these image sets turn out to be closed and convex subsets of 
^NxN- hi fact we claim that the following set is the image set of 4>: 

^NxN-= {Ae^£ NxN ;A T =A, A>0, A1 = 0} . 

To begin with it is clear that the image of 4> is equal to the image of $ and that it 
is included in &nxN, i-e. 

(RP x---RP)/E(p) 3 ([x 1 ,x 2 ,...,x N ])^C(xi,x 2 ,...,x N ) G &> NxN . 

The following proposition implies that £?NxN is the image set of 4> and it is the 
main result of this subsection. 
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Proposition 2 The mapping <J>: 

(R p x---xRP)/E(p) 3 ([xi,x 2 ,...,x N ]) M-C(xi,x 2 ,...,x N ) G &> NxN 
is a bijection. 

Proof. If A G &NxN we can > by the spectral theorem, find a unique symmet- 
ric and positive N x N matrix B = [bi,D2, . . . ,bjv] (the square root of A) with 
rank(B) = rank(A) such that B 2 = A and B 1 = 0. 

We now map the points (bi,D2, • • • ,b#) laying in R w isometrically to points 
(xi,X2, . . . ,xjy) in R p . This is trivially possible since p > N. The corresponding 
covariance matrix C(xi, . . . ,x#) will be equal to A. This proves surjectivity. That 
the mapping 4> is injective follows directly from the following lemma. 

Lemma 2.2 Let {yfc}f =1 and {y^}^ =1 be two sets of vectors in W. If 

ylyj = fkfj f°r j,k= 1,2,..., N, 

then there exists an S G O(p) such that 

S{y k )=yk fork=\,2,...,N. 

Proof. Use the Gram-Schmidt orfhogonalization procedure on both sets at the 
same time, q.e.d. 
q.e.d. 

We will in our explorations of high dimensional real world datasets below 
use "artificial" distance matrices constructed from geodesic distances on carefully 
created graphs connecting the samples or the variables. These distance matrices 
are converted to unique corresponding covariance matrices which in turn, as de- 
scribed above, give rise to canonical elements in (R p x ■■ ■ x ~R P )/E(p). We then 
pick sample centered representatives on which we perform PCA. In this way we 
can visualize low dimensional "approximative graph distances" in the dataset. Us- 
ing graphs in the sample set constructed from a k nearest neighbors or a locally 
euclidean approximation procedure, this approach corresponds to the ISOMAP 
algorithm introduced by Tenenbaum et. al. [52]. The ISOMAP algorithm can as 
we will see below be very useful in the exploration of DNA microarray data, see 
Nilsson et. al. Il38l for one of the first applications of ISOMAP in this field. 

We finally remark that if a proposed artificial distance or covariance matrix 
does not have the correct structure, i.e. if for example a proposed covariance 
matrix does not belong to &NxN, we begin by projecting the proposed covariance 
matrix onto the unique nearest point in the closed and convex set ^MxN an d then 
apply the scheme presented above to that point. 
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3 The basic statistical framework. 



We will here fix some notation and for the non- statistician reader's convenience at 
the same time recapitulate some standard multivariate statistical theory. In partic- 
ular we want to stress some basic facts concerning robustness of statistical testing. 

Let y be the sample space consisting of all possible samples (in our example 
datasets equal to all trials of patients) equipped with a probability measure P : 
2^ — y [0, +°°] and let X = (X\,. . . ,X p ) T be a random vector from y into R p . 
The coordinate functions X ; : y — y R, i = 1 , 2, . . . , p are random variables and in 
our example datasets they represent the expression levels of the different genes. 

We will be interested in the law of X, i.e. the induced probability measure 
P(X~ l (•)) defined on the measurable subsets of R p . If it is absolutely continuous 
with respect to Lebesgue measure then there exists a probability density function 
(pdf) fx(-) : R p — ► [0,oo) that belongs to Jz? l (R p ) and satisfies 

P({sey;X{s)eA})=jj x (x)dx (3.1) 

for all events (i.e. all Lebesgue measurable subsets) A C R p . This means that the 
pdf fx(-) contains all necessary information in order to compute the probability 
that an event has occurred, i.e. that the values of X belong to a certain given set 
A c R p . 

All statistical inference procedures are concerned with trying to learn as much 
as possible about an at least partly unknown induced probability measure P(X~ l (■)) 
from a given set of N observations, {x^x 2 , . . . ,x N } (withx* 6 R p for i = 1,2, . . . , AT), 
of the underlying random vector X. 

Often we then assume that we know something about the structure of the cor- 
responding pdf fx(-) and we try to make statistical inferences about the detailed 
form of the function fx(-). 

The most important probability distribution in multivariate statistics is the 
multivariate normal distribution. In R p it is given by the p-variate pdf n : R p — y 
(0,°o) where 

n{x) := (2^)-W2| r |-V2 e -i(x-M) r r- 1 (^) ;xG rp. q.2) 

It is characterized by the symmetric and positive definite p x p matrix T and the 
p-column vector ji, and |T| stands for the absolute value of the determinant of T. 



If a random vector X : y — y R p has the p-variate normal pdf (3.2 ) we say that X 
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has the N(jU,T) distribution. If X has the N(/i,T) distribution then the expected 
value of X is equal to /i i.e. 

S{X) := I X(s)dP = ji, (3.3) 

and the covariance matrix of X is equal to T, i.e. 

<g(X):= I \x-S{X)){X-S{X)) J dP = T. (3.4) 
Js 

Assume now thatX 1 ,* 2 ,...,^ are given independent and identically distributed 
(i.i.d.) random vectors. A test statistic J? is then a function (X l ,X 2 , . . . ,X N ) 
2?(X l ,X 2 , . . . ,X N ). Two important test statistics are the sample mean vector of 
a sample of size N 



1 '» 

M I— i ' 



and the sample covariance matrix of a sample of size Af 

1 N 



— ^(r-x^r-xy 



IfX\X 2 ,...,X N are independent and N(jU,T) distributed, then the mean X 
has the N(/l,^r) distribution. In fact this result is asymptotically robust with 
respect to the underlying distribution. This is a consequence of the well known 
and celebrated central limit theorem: 

Theorem 3.1 If the random p vectors X { ,X 2 ,X 3 , . . . are independent and identi- 
cally distributed with means /i G R p and covariance matrices T, then the limiting 
distribution of 

(AO 1/2 (^-M) 

asN — >-oowN(0,r). 

The central limit theorem tells us that, if we know nothing and still need to assume 
some structure on the underlying p.d.f., then asymptotically the N(fl,T) distribu- 
tion is the only reasonable assumption. The distributions of different statistics are 
of course more or less sensitive to the underlying distribution. In particular the 
standard univariate Student t-statistic, used to draw inferences about a univariate 
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sample mean, is very robust with respect to the underlying probability distribu- 
tion. In for example the study on statistical robustness pT| the authors conclude 
that: 

"...the two-sample t-test is so robust that it can be recommended in nearly all 
applications." 

This is in contrast with many statistics connected with the sample covariance 
matrix. A central example in multivariate analysis is the set of eigenvalues of 
the sample covariance matrix. These statistics have a more complicated behavior. 
First of all, if X\X 2 , ...,X N with values in are independent and N(/i,r) dis- 
tributed then the sample covariance matrix is said to have a Wishart distribution 
W p (N 1 T).lfN>p the Wishart distribution is absolutely continuous with respect 
to Lebesgue measure and the probability density function is explicitly known, see 
e.g. Theorem 7.2.2. in J3]|. If N >> p then the eigenvalues of the sample covari- 
ance matrix are good estimators for the corresponding eigenvalues of the under- 
lying covariance matrix T, see [0 and 0. In the applications we have in mind 
we often have the reverse situation, i.e. p >> TV, and then the eigenvalues for the 
sample covariance matrix are far from consistent estimators for the corresponding 
eigenvalues of the underlying covariance matrix. In fact if the underlying covari- 
ance matrix is the identity matrix it is known (under certain growth conditions on 
the underlying distribution) that if we let p depend on /V and if p/N — > y G (0, °°) 
as N — > oo, then the largest eigenvalue for the sample covariance matrix tends to 
(1 + tJy) 2 , see e.g. [55 1, and not to 1 as one maybe could have expected. This 
result is interesting and can be useful, but there are many open questions concern- 
ing the asymptotic theory for the "large p, large N case", in particular if we go 
beyond the case of normally distributed data, see e.g. 0, lfT8ll . 11261 . [|27l and 
113011 for an overview of the current state of the art. To estimate the information 
content or signal to noise ratio in our PCA plots we will therefore rely mainly on 
randomization tests and not on the (not well enough developed) asymptotic theory 
for the distributions of eigenvalues of random matrices. 

4 Controlling the false discovery rate 

When we perform for example a Student t-test to estimate whether or not two 
groups of samples have the same mean value for a specific variable we are per- 
forming a hypothesis test. When we do the same thing for a large number of 
variables at the same time we are testing one hypothesis for each and every vari- 
able. It is often the case in the applications we have in mind that tens of thousands 
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of features are tested at the same time against some null hypothesis Hq, e.g. that 
the mean values in two given groups are identical. To account for this multi- 
ple hypotheses testing, several methods have been proposed, see e.g. Il54l for an 
overview and comparison of some existing methods. We will give a brief review 
of some basic notions. 

Following the seminal paper by Benjamini and Hochberg lfT2ll . we introduce 
the following notation. We consider the problem of testing m null hypotheses Hq 
against the alternative hypothesis H\. We let niQ denote the number of true nulls. 
We then let R denote the total number of rejections, which we will call the total 
number of statistical discoveries, and let V denote the number of false rejections. 
In addition we introduce stochastic variables U and T according to Table [T] 

Table 1: Test statistics 





Accept Hq 


Reject Hq 


Total 


Hq true 


U 


V 


m 


H\ true 


T 


S 


mi 




m — R 


R 


m 



The false discovery rate was loosely defined by Benjamini and Hochberg as 
the expected value E{^). More precisely the false discovery rate is defined as 

FDR:=E(-\R>0)P(R>0). (4.1) 
R 

The false discovery rate measures the proportion of Type I errors among the sta- 
tistical discoveries. Analogously we define corresponding statistics according to 
Table [2} We note that the FNDR is precisely the proportion of Type II errors 

Table 2: Statistical discovery rates 

Expected value Name 

E(V/R) False discovery rate (FDR) 

E(T/(m — R)) False negative discovery rate (FNDR) 
E(T/(T + S)) False negative rate (FNR) 

E(V/(U + V)) False positive rate (FPR) 



among the accepted null hypotheses , i.e. the non-discoveries. In the datasets that 
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we encounter within bioinformatics we often suspect mi «m and so if R, which 
we can observe, is relatively small, then the FNDR is controlled at a low level. As 
pointed out in [|39ll , apart from the FDR which measures the proportion of false 
positive discoveries, we usually are interested in also controlling the FNR, i.e. we 
do not want to miss too many true statistical discoveries. We will address this 
by using visualization and knowledge based evaluation to support the statistical 
analysis. 

In our exploration scheme presented in the next section we will use the step 
down procedure on the entire list of p-values for the statistical test under consid- 
eration suggested by Benjamini and Hochberg in [fT2l to control the FDR. We will 
also use the q- value, computable for each separate variable, introduced by Storey, 
see [48 J and [|49l . The q- value in our analyses is defined as the lowest FDR for 
which the particular hypothesis under consideration would be accepted under the 
Benjami-Hochberg step down procedure. A practical and reasonable threshold 
level for the q- value to be informative that we will use is q < 0.2. 

5 The basic exploration scheme 

We will look for significant signals in our data set in order to use them e.g. as a 
basis for variable selection, sample classification and clustering. 

If there are enough samples one should first of all randomly partition the set of 
samples into a training set and a testing set, perform the analysis with the training 
set and then validate findings using the testing set. This should then ideally be 
repeated several times with different partitions. With very few samples this is 
not always feasible and then, in addition to statistical measures, one is left with 
using knowledge based evaluation. One should then remember that the ultimate 
goal of the entire exploration is to add pieces of new knowledge to an already 
existing knowledge structure. To facilitate knowledge based evaluation, the entire 
exploration scheme is throughout guided by visualization using PCA biplots. 

When looking for significant signals in the data, one overall rule that we follow 

is: 

• Detect and then remove the strongest present signal. 

Detect a signal can e.g. mean to find a sample cluster and a connected list of 
variables that discriminate the cluster. We can then for example (re)classify the 
sample cluster in order to use the new classification to perform more statistical 
tests. After some statistical validation we then often remove a detected signal, 
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e.g. a sample cluster, in order to avoid that a strong signal obscures a weaker but 
still detectable signal in the data. Sometimes it is of course convenient to add the 
strong signal again at a later stage in order to use it as a reference. 

We must constantly be aware of the possibility of outliers or artifacts in our 
data and so we must: 

• Detect and remove possible artifacts or outliers. 

An artifact is by definition a detectable signal that is unrelated to the basic mech- 
anisms that we are exploring. An artifact can e.g. be created by different experi- 
mental setups, resulting in a signal in the data that represents different experimen- 
tal conditions. Normally if we detect a suspected artifact we want to, as far as 
possible, eliminate the influence of the suspected artifact on our data. When we 
do this we must be aware that we normally reduce the degrees of freedom in our 
data. The most common case is to eliminate a single nominal factor resulting in 
a splitting of our data in subgroups. In this case we will mean-center each group 
discriminated by the nominal factor, and then analyze the data as usual, with an 
adjusted number of degrees of freedom. 

The following basic exploration scheme is used 

• Reduce noise by PCA and variance filtering. Assess the signal/noise ratio 
in various low dimensional PCA projections and estimate the projection 
information contents by randomization. 

• Perform statistical tests. Evaluate the statistical tests using the FDR, ran- 
domization and permutation tests. 

• Use graph-based multidimensional scaling (ISOMAP) to search for sig- 
nals/clusters. 

The above scheme is iterated until "all" significant signals are found and it is 
guided and coordinated by synchronized PCA-biplot visualizations. 

6 Some biological background concerning the ex- 
ample datasets. 

The rapid development of new biological measurement methods makes it possi- 
ble to explore several types of genetic alterations in a high-throughput manner. 
Different types of microarrays enable researchers to simultaneously monitor the 
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expression levels of tens of thousands of genes. The available information content 
concerning genes, gene products and regulatory pathways is accordingly grow- 
ing steadily. Useful bioinformatics databases today include the Gene Ontology 
project (GO) [5] and the Kyoto Encyclopedia of Genes and Genomes (KEGG) 
[|28ll which are initiatives with the aim of standardizing the representation of genes, 
gene products and pathways across species and databases. A substantial collec- 
tion of functionally related gene sets can also be found at the Broad Institute's 
Molecular Signatures Database (MSigDB) [36J together with the implemented 
computational method Gene Set Enrichment Analysis (GSEA) ll37l . ll5Tll . The 
method GSEA is designed to determine whether an a priori defined set of genes 
shows statistically significant and concordant differences between two biological 
states in a given dataset. 

Bioinformatic data sets are often uploaded by researchers to sites such as the 
National Centre for Biotechnology Information's database Gene Expression Om- 
nibus (GEO) [fT9l or to the European Bioinformatics Institute's database Array- 
Express BH. In addition data are often made available at local sites maintained by 
separate institutes or universities. 

7 Analysis of microarray data sets 

There are many bioinformatic and statistical challenges that remain unsolved or 
are only partly solved concerning microarray data. As explained in [|43l , these in- 
clude normalization, variable selection, classification and clustering. This state of 
affairs is partly due to the fact that we know very little in general about underlying 
statistical distributions. This makes statistical robustness a key issue concerning 
all proposed statistical methods in this field and at the same time shows that new 
methodologies must always be evaluated using a knowledge based approach and 
supported by accompanying new biological findings. We will not comment on the 
important problems of normalization in what follows but refer to e.g. (6) where 
different normalization procedures for the Affymetrix platforms are compared. In 
addition, microarray data often have a non negligible amount of missing values. 
In our example data sets we will, when needed, impute missing values using the 
K- nearest neighbors method as described in I153L All visualizations and analyses 
are performed using the software Qlucore Omics Explorer ll45l . 
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7.1 Effects of cigarette smoke on the human epithelial cell tran- 
scriptome. 

We begin by looking at a gene expression dataset coming from the study by Spira 
et. al. Il46l of effects of cigarette smoke on the human epithelial cell transcrip- 
tome. It can be downloaded from National Center for Biotechnology Informa- 
tions (NCBI) Gene Expression Omnibus (GEO) (DataSet GDS534, accession no. 
GSE994). It contains measurements from 75 subjects consisting of 34 current 
smokers, 18 former smokers and 23 healthy never smokers. The platform used to 
collect the data was Affymetrix HG-U133A Chip using the Affymetrix Microar- 
ray suite to select, prepare and normalize the data, see 11461 for details. 

One of the primary goals of the investigation in [|46ll was to find genes that are 
responsible for distinguishing between current smokers and never smokers and 
also investigate how these genes behaved when a subject quit smoking by looking 
at the expression levels for these genes in the group of former smokers. We will 
here focus on finding genes that discriminate the groups of current smokers and 
never smokers. 

• We begin our exploration scheme by estimating the signal/noise ratio in 
a sample PCA projection based on the three first principal components. 

We use an SVD on the data correlation matrix, i.e. the covariance matrix for the 
variance normalized variables. In Figure [T] we see the first three principal com- 
ponents for ImL and the 75 patients plotted. The first three principal components 
contain 25% of the total variance in the dataset and so for this 3-D projection 
«2({ 1,2, 3}, obsr) = 0.25. Using randomization we estimate the expected value 
for a corresponding dataset (i.e. a dataset containing the same number of samples 
and variables) built on independent and normally distributed variables to be ap- 
proximately 0C2({l 5 2,3},ran<i) = 0.035. We have thus captured around 7 times 
more variation than what we would have expected if the variables were inde- 
pendent and normally distributed. This indicates that we do have strong signals 
present in the dataset. 

• Following our exploration scheme we now look for possible outliers and 
artifacts. 

The projected subjects are colored according to smoking history, but it is clear 
from Figure[T]that most of the variance in the first principal component (containing 
13% of the total variance in the data) comes from a signal that has a very weak 
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Figure 1: The 34 current smokers (red), 18 former smokers (blue) and 23 never 
smokers (green) projected onto the three first principal components. The separa- 
tion into two groups is not associated with any supplied clinical annotation and is 
thus a suspected artifact. 




Figure 2: The red samples have high description numbers (> 58) and the green 
samples have low description numbers (< 54). The blue sample has number 5. 

association with smoking history. We nevertheless see a clear splitting into two 
subgroups. Looking at supplied clinical annotations one can conclude that the two 
groups are not associated to gender, age or race traits. Instead one finds that all 
the subjects in one of the groups have low subject description numbers whereas all 
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the subjects except one in the other group have high subject description numbers. 
In Figure [2] we have colored the subjects according to description number. This 
suspected artifact signal does not correspond to any in the dataset (Dataset GDS 
534, NCBIs GEO) supplied clinical annotation. One can hypothesize that the 
description number could reflect for instance the order in which the samples were 
gathered and thus could be an artifact. Even if the two groups actually correspond 
to some interesting clinical variable, like disease state, that we should investigate 
separately, we will consider the splitting to be an artifact in our investigation. We 
are interested in using all the assembled data to look for genes that discriminate 
between current smokers and never smokers. We thus eliminate the suspected 
artifact by mean-centering the two main (artifact) groups. After elimination of the 
strong artifact signal, the first three principal components contain "only" 17% of 
the total variation. 

• Following our exploration scheme we filter the genes with respect to 
variance visually searching for a possibly informative three dimensional 
projection. 




Figure 3: We have filtered by variance keeping the 630 most variable genes. It 
is interesting to see that the third principle component containing 9% of the total 
variance separates the current smokers (red) from the never smokers (green) quite 
well. 

When we filter down to the 630 most variable genes, the three first princi- 
pal components have an L 2 -projection content of 0t2({l,2,3}) = 0.42, whereas 
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Table 3: Top genes upregulated in the current smokers group and downregulated 
in the never smokers group. 



Gene symbol 


g-value 


NOD1 


5 59 1 0406777 1 R94e-0R 


GPX9 


9 11 149939391 979e-07 


AT DH3A1 


9 11 149939391 979e-07 


fT DN10 


1 45691 4391 69953e-06 


FTH1 


4.72936617815058e-06 


TALDOl 


4.72936617815058e-06 


TXN 


4.72936617815058e-06 


MUC5AC 


3.77806345774405e-05 


TSPAN1 


4.50425200297664e-05 


PRDX1 


4.58227420582093e-05 


MUC5AC 


4.99131989472012e-05 


AKR1C2 


5.72678146958168e-05 


CEACAM6 


0.000107637125805187 


AKR1C1 


0.000195523829628407 


TSPAN8 


0.000206106293159401 


AKR1C3 


0.000265342898771159 



by estimation using randomization we would have expected it to be 0.065. The 
projection in Figure [3] is thus probably informative. We have again colored the 
samples according to smoking history as above. The third principal component, 
containing 9% of the total variance, can be seen to quite decently separate the cur- 
rent smokers from the never smokers. We note that this was impossible to achieve 
without removing the artifact signal since the artifact signal completely obscured 
this separation. 

• Using the variance filtered list of 630 genes as a basis, following our 
exploration scheme, we now perform a series of Student t-tests between 
the groups of current smokers and never smokers, i.e. 34 + 23 = 57 
different subjects. 

For a specific level of significance we compute the 3-dimensional (i.e. we let 
S = {1,2,3}) L 2 -projection content resulting when we keep all the rejected null 
hypotheses, i.e. statistical discoveries. For a sequence of t-tests parameterized by 
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Table 4: Top genes downregulated in the current smokers group and upregulated 
in the never smokers group. 



Gene symbol 


<^-value 


MT1G 


4.03809377378893e-07 


MT1X 


4.72936617815058e-06 


MUC5B 


2.38 1989034023 17e-05 


CD81 


3.1 60522 1864278e-05 


MT1L 


3.1 60522 1864278e-05 


MT1H 


3.1 60522 1864278e-05 


SCGB1A1 


4.50425200297664e-05 


EPAS1 


4.63861480935914e-05 


FABP6 


0.00017793865432854 


MT2A 


0.000236481909692626 


MT1P2 


0.000251264650053933 




Figure 4: A synchronized biplot showing samples to the left and variables to the 
right. 



the level of significance we now try to find a small level of significance and at 
the same time an observed L 2 -projection content with a large quotient compared 
to the expected projection content estimated by randomization. We supervise this 
procedure visually using three dimensional PCA-projections looking for visually 
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clear patterns. For a level of significance of 0.00005, leaving a total of 43 genes 
(rejected nulls) and an FDR of 0.0007 we have (^({l, 2, 3},obsr) = 0.71 whereas 
the expected projection content for randomized data 0C2({1,2, 3}, rand) =0.21. 
We have thus captured more than 3 times of the expected projection content and 
at the same time approximately 0.0007 x 43 = 0.0301 genes are false discoveries 
and so with high probability we have found 43 potentially important biomarkers. 
We now visualize all 75 subjects using these 43 genes as variables. In Figure|4]we 
see a synchronized biplot with samples to the left and variables to the right. The 
sample plot shows a perfect separation of current smokers and never smokers. In 
the variable plot we see genes (green) that are upregulated in the current smokers 
group to the far right. The top genes according to g-value for the Student t-test 
between current smokers and never smokers, that are upregulated in the current 
smokers group and downregulated in the never smokers group, are given in Table 
[3] In Table [4] we list the top genes that are downregulated in the current smokers 
group and upregulated among the never smokers. 

7.2 Analysis of various muscle diseases. 

In the study by Bakay et. al. [ITOl the authors studied 125 human muscle biop- 
sies from 13 diagnostic groups suffering from various muscle diseases. The plat- 
forms used were Affymetrix U133A and U133B chips. The dataset can be down- 
loaded from NCBIs Gene Expression Omnibus (DataSet GDS2855, accession no. 
GSE3307). We will analyze the dataset looking for phenotypic classifications and 
also looking for biomarkers for the different phenotypes. 

• We first use 3 -dimensional PCA-projections of the samples of the data 
correlation matrix, filtering the genes with respect to variance and vi- 
sually searching for clear patterns. 

When filtering out the 300 genes having most variability over the sample set 
we see several samples clearly distinguishing themselves and we capture 46% 
of the total variance compared to the, by randomization estimated, expected 6%. 
The plot in Figure|5]thus contains strong signals. Comparing with the color legend 
we conclude that the patients suffering from spastic paraplegia (Spg) contribute a 
strong signal. More precisely, three of the subjects suffering from the variant Spg- 
4 clearly distinguish themselves, while the remaining patient in the Spg-group 
suffering from Spg-7 falls close to the rest of the samples. 

• We perform Student t-tests between the spastic paraplegia group and 
the normal group. 
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I I Becker muscular dystrophy 
I I Calpain 3 mutation 
| Duchenne muscular dystrophy 
| Emery -Dreifuss muscular dystrophy 
| EmeryDreifussFSHD 
□ FKRP mutation 
HFSH Dysferlin 

Hfsh-dmd 

| acute quadriplegic myopathy 
1 amyotophic lateral sclerosis 
I I dysferlin mutation 
CD juvenile dermatomyositis 
I I normal 

I I spastic paraplegia 

(a) A three dimensional PCA plot capturing 46% (b) Color legend 

of the total variance. 

Figure 5: PCA-projection of samples based on the variance filtered top 300 genes. 
Three out of four subjects in the group Spastic paraplegia (Spg) clearly distinguish 
themselves. These three suffer from Spg-4, while the remaining Spg-patient suf- 
fers from Spg-7. 
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Table 5: Top genes upregulated in the Spastic paraplegia group (Spg-4). 



Gene symbol 


g-value 


RAB40C 


0.0000496417 


SFXN5 


0.000766873 


CLPTM1L 


0.00144164 


FEM1A 


0.0018485 


HDGF2 


0.00188435 


WDR24 


0.00188435 


NAPSB 


0.00188435 


ANKRD23 


0.00188435 



As before we now, in three dimensional PCA-projections, visually search for 
clearly distinguishable patterns in a sequence of Student t-tests parametrized by 
level of significance, while at the same time trying to obtain a small FDR. At a 
level of significance of 0.00001, leaving a total of 37 genes (rejected nulls) with 
an FDR of 0.006, the first three principal components capture 81% of the variance 
compared to the, by randomization, expected 31%. Table [5] lists the top genes 
upregulated in the group spastic paraplegia. We can add that these genes are all 
strongly upregulated for the three particular subjects suffering from Spg-4, while 
that pattern is less clear for the patient suffering from Spg-7. 

• In order to find possibly obscured signals, we now remove the Spastic 
paraplegia group from the analysis. 

We also remove the Normal group from the analysis since we really want to com- 
pare the different diseases. Starting anew with the entire set of genes, filtering 
with respect to variance, we visually obtain clear patterns for the 442 most vari- 
able genes. The first three principal components capture 46% of the total variance 
compared to the, by randomization estimated, expected 6%. 

• Using these 442 most variable genes as a basis for the analysis, we now 
construct a graph connecting every sample with its two nearest (using 
euclidean distances in the 442-dimensional space) neighbors. 

As described in the section on multidimensional scaling above, we now com- 
pute geodesic distances in the graph between samples, and construct a resulting 
distance (between samples) matrix. We then convert this distance matrix to a 
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corresponding covariance matrix and finally perform a PCA on this covariance 
matrix. The resulting plot (together with the used graph) of the so constructed 
three dimensional PCA-projection is depicted in Figure [6j 




Figure 6: Effect of the ISOMAP-algorithm. We can identify a couple of clusters 
corresponding to the groups juvenile dermatomyositis, amyotophic lateral sclero- 
sis, acute quadriplegic myopathy and Emery Dreifuss FSHD. 

Comparing with the Color legend in Figure [5j we clearly see that the groups 
juvenile dermatomyositis, amyotophic lateral sclerosis, acute quadriplegic my- 
opathy and also Emery-Dreifuss FSHD distinguish themselves. 

One should now go on using Student t-tests to find biomarkers (i.e. genes) dis- 
tinguishing these different groups of patients, then eliminate these distinct groups 
and go on searching for more structure in the dataset. 

7.3 Pediatric Acute Lymphoblastic Leukemia (ALL) 

We will finally analyze a dataset consisting of gene expression profiles from 132 
different patients, all suffering from some type of pediatric acute lymphoblastic 
leukemia (ALL). For each patient the expression levels of 22282 genes are ana- 
lyzed. The dataset comes from the study by Ross et. al. Il44l and the primary 
data are available at the St. Jude Children's Research Hospital's website ll50l . 
The platform used to collect this example data set was Affymetrix HG-U133 chip, 
using the Affymetrix Microarray suite to select, prepare and normalize the data. 
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As before we start by performing an SVD on the data correlation matrix vi- 
sually searching for interesting patterns and assessing the signal to noise ratio by 
comparing the actual L 2 -projection content in the real world data projection with 
the expected L 2 -projection content in corresponding randomized data. 

• We filter the genes with respect to variance, looking for strong signals. 




(a) Projection capturing 38% of the total vari- (b) Color legend 
ance. 



Figure 7: Variance filtered PCA-projection of the correlation datamatrix based on 
873 genes. The group T-ALL clearly distinguish itself. 

In Figure [7] we see a plot of a three dimensional projection using the 873 
most variable genes as a basis for the analysis. We clearly see that the group T- 
ALL is mainly responsible for the signal resulting in the first principal component 
occupying 18% of the total variance. In fact by looking at supplied annotations 
we can conclude that all of the other subjects in the dataset are suffering from 
B-ALL, the other main ALL type. 

• We now perform Student t-tests between the group T-ALL and the rest. 
We parametrize by level of significance and visually search for clear 
patterns. 

In Figure [8] we see a biplot based on the 70 genes that best discriminate be- 
tween T-ALL and the rest. The FDR is extremely low FDR = 1. 13e — 24 telling 
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Figure 8: FDR = 1.13e — 24. A synchronized biplot showing samples to the left 
and genes to the right. The genes are colored according to their expression level 
in the T-ALL group. Red = upregulated and green = downregulated. 

us that with a very high probability the genes found are relevant discoveries. The 
most significantly upregulated genes in the T-ALL group are 

CD3D, CD7, TRD@, CD3E, SH2D1A and TRA@. 
The most significantly downregulated genes in the T-ALL group are 

CD74, HLA-DRA, HLA-DRB, HLA-DQB and BLNK. 

By comparing with gene-lists from the MSig Data Base (see I; 361 ) we can 
see that the genes that are upregulated in the T-ALL group (CD3D, CD7 and 
CD3E) are represented in lists of genes connected to lymphocyte activation and 
lymphocyte differentiation. 

• We now remove the group T-ALL from the analysis and search for visu- 
ally clear patterns among three dimensional PCA-projections filtrating 
the genes with respect to variance. 

Starting anew with the entire list of genes, filtering with respect to variance, a clear 
pattern is obtained for the 226 most variable genes. We capture 43% of the total 
variance as compared to the expected 6.5%. We thus have strong signals present. 

• Using these 226 most variable genes as a basis for the analysis, we now 
construct a graph connecting every sample with its two nearest neigh- 
bors. 
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Figure 9: Effect of the ISOMAP-algorithm. We can clearly see the groups E2A- 
PBX1, MLL and TEL-AML1. The group TEL-AML1 is connected to a subgroup 
of the group called Other (white). 

We now perform the ISOMAP-algorithm with respect to this graph. The resulting 
plot (together with the used graph) of the so constructed three dimensional PCA- 
projection is depicted in Figure |9j We can clearly distinguish the groups E2A- 
PBX1, MLL and TEL-AML1. The group TEL-AML1 is connected to a subgroup 
of the group called Other. This subgroup actually corresponds to the Novel Group 
discovered in the study by Ross et.al. 11441 . Note that by using ISOMAP we 
discovered this Novel subgroup only by variance filtering the genes showing that 
ISOMAP is a useful tool for visually supervised clustering. 
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